NOTICE 


THIS DOCUMENT HAS BEEN REPRODUCED FROM 
MICROFICHE. ALTHOUGH IT IS RECOGNIZED THAT 
CERTAIN PORTIONS ARE ILLEGIBLE, IT IS BEING RELEASED 
IN THE INTEREST OF MAKING AVAILABLE AS MUCH 
INFORMATION AS POSSIBLE 






•Made available under HA:-' 
in the interest ol early and «r 
semination ol Earth Resources 
Prograff information and v:itno .» 
for anv use made thereof 














TELL. US-NEWSLETTER 11 


November 1979 


a 0 - 10333 


FURTHER DEVELOPMENTS OF THE TELL-US MODEL : 


I AN IMPLICIT FINITE DIFFERENCE SCHEME FOR THE NUMERICAL 
APPROXIMATION OF THE GROUND HEAT FLUX 

II A SIMPLE ALGORITHM FOR ESTIMATING THE ACTUAL AND 
POTENTIAL EVAPOTRANSP I RATION OF VEGETATED SURFACES 
FROM ONE REMOTELY SENSED SURFACE TEMPERATURE NEAR 
THE DAILY MAXIMUM 


by : J. HUY GEN 

Joint Research Centre, Ispra Establishment 


Correspondence regarding this series of newsletters should be addressed to : 

P. REIN/GER Commission of the European Communities 
JOINT RESEARCH CENTRE 
1-21020 Ispra (Varese) - Italy 


AN IMPLICIT FINITE DIFFERENCE SCHEME FOR THE 
NUMERICAL APPROXIMATION OF THE GROUND HEAT FLUX 


INTRODUCTION 

In the appendix of TELLUS-Nowsletter 8, May 1979^) it had been 
shown that using an explicit finite difference scheme to approximate 
the flow equation for heat transport in soil, the time- and depth step 
must be kept small as the stability criterion may not be exceeded 
and the approximation may not differ too much from the exact solu- 
tion. It also appeared that the scheme of Du Fort-Frankel is only 
slightly superior in both respects to the classical explicit difference 
method. 


Implicit finite difference methods are known to be superior to explicit 
methods, especially with respect to stability. Therefore, they permit 
large time steps, keeping computer calculation time low. Good results 
have recently been obtained by building an implicit system into the 
TELL -US algorithm. 

THE IMPLICIT FINITE DIFF E RKNCE METHOD APPLIED 


To the one -dimensional he at flow equation (J) at Lime l and depth 
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the same equation at time t + ^t must be added to obtain: 
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According to C ranck and Nicolsm. cm left-hand side of i£q. (2) 
may be approximated by: 


A t 


(T(x,t + At) - T(x, t)), 


neglecting the derivatives of order two and higher. 

Writing also the right-hand side of Eq. (2) in difference form, the 
following scheme can be obtained: 
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— (T(x, t+At) - T(x, t)) * [t(x-Ax, t) - 2T(x,t) + T(x+Ax, t) ] 

A Ax 
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For .i grid expanding in the x-direetion Kq, (i) can be written »u<: 
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A i • time step (sec) 

A x depth atop, sjxieing of the upper two grid point a (m) 
KXF - expansion factor 

P (2*1-1), Id for the first grid point below the surface. 


To solve Kq. (4) u set of equations is needed. To give an example, 
let us assume a soil system divided into four layers, i. e. there are 
(4 tl) grid points from the surface (x 0) to the lower boundary depth 
(x d). ' 

When Kq. (4) is written out for successively x 1, l a set of (4-1) 
equations is formed, which reads in matrix notation; 
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Vector T(x, t) is the known temperature profile at time t, 

T(4, t + A t) * T(4, l) is the lower boundary temperature, 

T(0, t i A t;) is at first approximated as; 2* T(0, t) - T(0,t - A t). 



Kvorytlung on the right -hand t> itlo of Kq, ( p) in known and results 
in n single vector, in liu> program called FRODd. Tho linear system 
that remains can easily bo solved by converting tho matrix in upper 
triangular 1'onii and back-substitution. Tho solution vector is the 
tempo rature profile in the soil at t !• A t . 

The temperature at x 0 and x ^ 1 are entered into the energy ha lance 
equation. The first one will be modified during the iterative process, 
While T(l,t ‘I At) is strongly influenced by the estimated T(U, t I At), it 
is clear that considerable errors may arise when the soil heat flux is 
calculated from the modified surface temperature and the fixed value 
at \ l. Therefore the estimated T(U, t 4 A t) will be replaced by the 
surface temperature that comes out of the iterative process and the 
whole procedure from Kq» ( f >) will he repeated once. 

RKK!!t/i\S 

>, Vi. t . s*^i:*** fc* •mwJtK, -» »*Sf 

In Tables l A and l it the results of calculations with the implicit ver- 
sion are shown. Tho data set obtained during the Joint Flight Kxperi- 
ment at tirendon I nderwoodj U. N, on September id-14, 1077 was used 
for these calculations. 

Table 1 A may he compared with tables 1 A and .’A of TKI.l.Utf News- 
letter 8, Resemblance is best with ruble dA of that paper. This means 
that there remain differences with what can he seen as the mosl pro 
oi.se calculations, i, e. Table 1 A of Newsletto r S. However, these dif- 
ferences are not dramatically largo. 

To recapitulate: a high heat capacity value of x l O' J/kg/K as ori- 
ginally applied in TKL1.-US is not acceptable as shown in Newsletter d, 
A realistic heal capacity value has to he accompanied by such a, small 
time step, in order to keep the explici . amorical seheate stable, that 
total Cvilculation time becomes ext rucruiuu rily high (alternative 1), 

With an explicit scheme and .(’“realistic l.eat capacity value some 
speeding up is possible by increas . v. c depth step to a round d cm 
and the time step to five m’.'. d • v ,.i. amative d)» The third alternative,- 

the proposed implicit sc name with a time step of one hour and a depth 
step of .1 cm, reduces calculation time with a factor four compared to 
alternative .1 (Id rosp, ‘40 see). 

It appeared from the treatment of several other data sets that TKl.l.-lUR 
becomes rat re r sensitive to heat capacity, when a realistic value is 
taken. Therefore it may he suggested to vary the heat capacity value 
with thermal inertia, which is from a physical point of view certainly 
defemlable. A look-up table calculated with a floating heat capacity 
value is shown as Table 111. 





TABLE 1A - Look-Up Table Calculated with the 
Version of the TELL-US >/ode 
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Programming Aspects of the Implicit Approach 
(The soil is divided into 10 layers between the surface and the 
r , lower boundary depth) 

C To calculate for every new thermal inertia value 

CAP a 2K6 (or a function of thermal inertia) 

COND - (THIN(J)±±2)/CAP 
Do 2 I - 1,9 

F0(I) e TS !bC ON D / C AP/( (E XF + 1 ) & DS ±*2 * EXF *±(2*1-1) 
KOI (il = -((F0(l)*EXF + F0(I) + l)/F0(I)/EXF) 

2 F0?.(l) « ((F0(I)±EXF + F0(I) - 1 )/F0(I)/EXF) 

C Description of matrix A (9*9) 

Do 11 1-1,9 
Do 11 J a 1, 9 
11 A(I, J)=0. 0 
Do 13 1=1,9 

13 A(I, I) = F01(I) 

Do .14 I = 1,8 

14 A(I, I + 1) a 1/EXF 
Do 1 5 I a 2, 9 

1 5 A(I, I - 1 ) = 1 

C Description of matrix B (9.4 9) 

Do 16 1=1, 9 
Do 17 J « 1,9 
17 B(I, J) S -A(I, J) 

16 B(I, 1) = F02(I) 

C --------------- 

C To calculate evex*y now time step: 

C matrix B multiplied with vector T(x, t) results in the vector 

C PRODl(x). 

C The vector T(x, t) is called TOLD x) 

Do 100 1*1,9 
PROD 1(1) --- 0. 0 
1 MIN 1 a 1-1 

IF (I MIN 1. LT. 1) I MIN 1 « 1 
I PLUS I =1+1 

IF(I PLUS 1. GT. 9) I PLUS 1 = 9 
Do 100 J a 1 MINI, IPLUS 1 
100 PROD 1(1) = PRODI(I) + B(I, J) ± TOLD(J) 


C calculation of PROD2(x) 

C T(0, t + A t) is called; STNEY7 

C T(0, t) is called: STOLD 

C T(0, t - M) is called; STVOLD 

CCCC STNEW is only a first approximation 

STNEW = 2 jfc STOLD - STVOLD 
PROD2(l) = PRODl(l) - STOLD - STNEW 
Do 200 I a 2, 8 
200 PROD 2(1) = PRODI(I) 

PROD2(9) = PROD 1(9) - 2 ifc LBT/EXF 

G calculation of new temperature profile 

C the solution vector T(x, t *1- At) is called TNEW(x) 

C COPY A IN D 

Do 300 1 = 1,9 
Do 300 J = 1,9 
300 D(I, J) = A(l, J) 

EPS = 1, Q E-10 
Do 306 I « 1,8 
Do 306 J = 1+1,9 

IF(ABS(D(J, I)). LT. EPS) GOTO 306 
HOLD = D( J, I) 

Do 307 K = 1, 9 

307 D(J, K) = D(I, K) - D(I, l)/HOLD* D(J, K) 

306 PROD2( J) = PROD 2 (I) - D(I, l)/HOLD S PROD2(J) 

TNEW(9) = PROD2(9)/D(9, 9) 

Do 308 IP = 1, 8 
I = 9-IP 
I PLUS 1=1+1 
SOM = 0. 0 

Do 309 J = I PLUS 1, 9 
309 SOM = SOM + D(l, J) ifc TNEW(J) 

308 TNEW(I) •- (PROD2(I) - SON;)/. j(x, .) 

All energy balance components now can be calculated as a function 
of first approximation STNEW. 

The modified STNEW coming out of the iterative process must be 
entered again in the equation: 

PROD2(l) = PROD 1(1) - STOLD - STNEW and the procedure to 
calculate vector TNEW must be repeated. 


A SIMPLE ALGORITHM FOR ESTIMATING THE ACTUAL AND 
POTENTIAL E V APOT RANSPIRAT ION OF VEGETATED SURFACES 
FROM ONE REMOTELY SENSED SURFACE TEMPERATURE NEAR 

THE DAILY MAXIMUM 


INTRODUCTION 


At present, and without considering other approaches or methods 
of calculation (SEGUIN, 197 8; JACKSON et al. , x 977), two models 
are available to the TELLUS project to calculate daily evapotran- 
s pi ration and soil moisture content from remotely sensed surface 
tern pe ratu res: 

1. The TERGIlA-model (SOER, 1977) was originally developed for 
grassland. It is still under investigation how the model could also 
be applied to other, more structured crops. 

1 . The TELL- US -model (ROSEMA and LYLE VELD, 1977) was set 
up for bare or scarcely vegetated soils, According to later de- 
velopments (KLAASSEN and ROSEMA, 1979), it should bo applicable 
in the same form also to vegetated surfaces. 

It has to be remarked that TERGRA in its basic form did not estimate 
soil moisture content (on the contrary: soil moisture potential at the 
beginning of simulation is an important input parameter), nor did it 
make use of a remotely sensed crop temperature. In later versions 
the model was adapted to the evaluation of soil moisture employing 
one remotely sensed surface temperature close to the daily maximum. 
For practical remote sensing purposes, the need for a large number 
of place -dependent input parameters still constitutes a drawback of 
the model. 

TELL-US uses two remotely sensed surface temperatures, close to 
the daily maximum and minimum temperatures, respectively. In this 
way, the course of the surface tec. v :. s rature through the day is more 
or less fixed. This method /.oreme to vary two parameters, surface 
relative humidity that predominantly determines the day temperature 
and thermal inertia that mainly influences the night temperature. 

The simulated daily course of the surface temperature is a function 
of these two variables. If, for a certain combination of these two 
variables, there is sufficient resemblance between the simulated 
and measured minimum and maximum surface temperatures, these 
values of surface relative humidity and thermal inertia are assumed 
to be valid for the system. Then, surface relative humidity permits 
the calculation of daily evapotranspiration axrd thermal inertia gives 
information about soil moisture content. 


i ^ 
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TELL- US needs only a few place - dependent input parameters and 
is therefore operationally more attractive than TERGRA. Its appli- 
cation to vegetated surfaces, as proposed by KLAASSKN and 
ROSEMA (1979), would therefore seem advantageous. Nevertheless, 
this application gives rise to a number of problems which will be 
discussed in the next section. A simplified algorithm is proposed 
which should avoid some of these problems. 

PROBLEMS RELATED TO THE APPLICATION OF T KLL-US TO 
VEGETATED SURFACES 


1 . During Daytime: 

The ground heat flux under a closed vegetative cover is known to 
be. small during daytime, about one order of magnitude less than 
net radiation. For its energy input, the soil depends on the vege- 
tation (radiative exchange between soil and crop) and on the rela- 
tively still air inside the canopy (conductive heat exchange). 

The soil surface will therefore usually be cooler than the vegeta- 
tive cover. In TELL-US, the cx-op temperature will be adapted as 
the upper boundax*y temperature of the soil system. This will 
result in an overestimate of the ground heat flux. As the surface 
temperature through the day is mor- or less fixed, an underesti- 
mate of the latent heat flux is inevitable, as net radiation and the 
sensible heat flux are fixed together with the surface temperature 
in the simulation process. 

1 . During Nighttime: 

At night, the vegetation also has an insulating effect and usually 
crop temperature will be lower than the soil surface temperature. 
Given the boundary condition i». TELL-US, again the soil heat flux 
is overestimated and, the crop temperature being fixed, the rmal 
inertia will be strongly underestimated. 

In the case of vegetated au./ faces, thermal inertia has exactly the 
same physical mc».-niu 0 as m the case of bare soils: it determines 
the thermal, properties of the soil and thus, indirectly, it gives 
information about soil moisture content. Contrary to bare soil, 
for a vegetated surface, also surface relative, humidity, or a 
correspondiixg parameter sxich as stomatal resistance, is corre- 
lated with bulk soil moisture content. It should be understood 
that in the case of bare soil, sxirface relative humidity and thermal 
inertia are independent variables while in the case of a vegetated 
surface they are strongly related. Therefore, it seems superfluous 
to distinguish between surface relative humidity and thermal inertia 
in the latter case and there remains little reason to use the nightly 
minimum surface temperature, as: 


\0 


aj the thermal inertia value obtained will give a wrong impression 
of soil moisture content, and 

b) wind speed is often low during the night which poses problems 
due to the great variability of the air temperature under these 
conditions (IIUYGEN and REINIGER, 1979). 


THE NEW APPROACH 

The simplified algorithm is a derivation of the TKLL-US model. 

The following simplifications have been introduced: 

1) Only the hours between sunrise and sunset are taken as the 
simulation period. 

<i) The ground heat flux is assumed to be 1 0 'jo of net radiation. 

3) Instead of surface relative humidity and thermal inertia, the course 
of the daily surface temperature is calculated for different but 
constant values of bulk stomatal resistance. 

For a certain stomatal resistance there will be agreement between 
the remotely sensed surface temperature close to the daily maxi- 
mum and the simulated crop temperature at that time. The stomatal 
resistance obtained in this way permits calculation of the daily eva- 
potrans pi ration. It will be shown that potential evapotrans piration 
can also be obtained by establishing a minimum value for the bulk 
stomatal resistance. 


The input parameters of the model are the. same as for the TELL-US 
model, i. e. daily maximum surface temperature, albedo, emissi- 
vity, crop height, slope dip and direction and hourly measurement 
of incoming short- and longwave rau...ion, dry- and wet bulb tempe- 
rature and windspeed. The incoming longwave radiation may even- 
tually be approximated by a Brum.-x.ype formula. 

A listing of the computer ..a is available on request. 

TEST OF THE ALGORITHM 


The algorithm was compared with results obtained with the TERGKA 
model for a simulated Super Test Data Set., given in Table 1. Re stilts 
of the test are shown in Fig. 1 with an example of the output of the 
model, given in Table 2. 

With the simplified algorithm potential evapotrans piration can be 
found by establishing a minimum value for bulk stomatal resistance. 
In the literature critical values for leaf w-ater potential can be found 
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up to which potential transpiration will be maintained. 

RIJTEMA and ABOUKHALED (1975) give a value of -1000 kPa for 
grass and wheat, whereas EHRLER et al. (1978) give a value of 
around -1500 kPa for wheat. 

With equation (1), used in TERGRA, leaf water potential may be 
translated into bulk stomatal resistance: 

RC -V 1/GIlC * (0. 05 * (PSI.L/100) ±± 2. 1 + 400/(RS + 1.5)) (1) 


where : 

RC = stomatal resistance (s/m) 

GHC = mean crop height (m) 

PSIL = leaf water potential (kPa) 

RS ~ incoming solar radiation (W/m^) 

For the Super Test Data Set, as shown in Fig. 1, the simplified 
algorithm calculates a total daily transpiration Of 4. 9 mm when the 
leaf water potential is kept at -1500 kPa. TERGRA also calculates 
a potential daily transpiration of 4. 9 mm (Fig. 1). For the data set 
in question one may conclude that the assumption of a critical leaf 
water potential of -1500 kPa leads to an estimate of potential evapo- 
transpiration comparable to that of the TERGRA model. However, it 
should be stressed that Eq. (1) was basically developed for grass. 

To be valid for crops with a different roughness, a suitable form of 
Eq. (1) should first be obtained. 

Actual evapotranspiration may be obtained by varying the bulk stomatal 
resistance until there is agreement between the simulated and the 
measured crop temperature. In Fig. 1 the actual evapotranspiration 
as a function of crop temperature at 13, 00 H, as calculated by 
TERGRA, is described by the full ctirve. The points marked X are 
values calculated by the simplified algorithm. For the atmospheric 
conditions of a rather high evapo muve demand, characteristic of 
this data set, the agree mem oetwoen the two models is rather good 
with a maximum difference of 5 °]o in the total daily evapotranspiration. 

CONCLUSION 

The proposed simplified evapotranspiration algorithm should be 
useful to estimate from remotely sensed surface temperatures both 
potential and actual evapotranspiration, important parameters in 
present day crop yield forecasting. It can, in principle, be applied 
to all crop types once the surface is completely covered. As with the 
TELL-US model a multidimensional look-up table can be constructed 


wluoh .should permit the mapping of evapotruuspi ration on .t pixel 
by pixel basis. 


Problems in the application of the algorithm are the estimation 
of roughness for different crop types and, probably, the an sump* 
lion of .1 bulk stomata! resistance constant during the simulation 
period. 

Likewise, the questions posed by the interpretation of remotely 
sensed surface temperatures of structured crops, still remain open. 


The algorithm 


should be further tested with experimental field data. 
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CROP TEMP. AT 1300 H (°C) 


Fig. 1 : Total daily evapotranspiration as a function of crop surface temperature 
at 13.00 H. 

Full curve : Calculated with TERGRA 

X • Calculated with the simplified model 




TABLE 1 - Super Test Data Set - Grassland 
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